Collision densities and mean residence times for d-dimensional exponential flights 
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In this paper we analyze some aspects of exponential flights, a stochastic process that governs the 
evolution of many random transport phenomena, such as neutron propagation, chemical/biological 
species migration, or electron motion. We introduce a general framework for d-dimensional setups, 
and emphasize that exponential flights represent a deceivingly simple system, where in most cases 
closed-form formulas can hardly be obtained. We derive a number of novel exact (where possible) or 
asymptotic results, among which the stationary probability density for 2d systems, a long-standing 
issue in Physics, and the mean residence time in a given volume. Bounded or unbounded, as well 
as scattering or absorbing domains are examined, and Monte Carlo simulations are performed so as 
to support our findings. 



I. INTRODUCTION 

Random walks are widely used in Physics so as to 
model the features of transport processes where the mi- 
grating (possibly massless) particle undergoes a series of 
random displacements as the effect of repeated collisions 
with the surrounding environment [l|-|4|. While much 
attention has been given to random walks on regular Eu- 
clidean lattices, and to the corresponding scaling limits, 
less has been comparatively devoted to the case where 
the direction of propagation can change continuously at 
each collision: for an historical review, see, e.g., [5|. Such 
processes, which are intimately connected to the Boltz- 
mann equation, have been named random flights, and 
play a prominent role in the description of, among oth- 
ers, neutron or photon propagation through matter [^-Q, 
chemical and biological species migration 9] , or electron 
motion in semiconductors (loj . 

Within the simplest formulation of this model, which 
was originally proposed by Pearson (1905) [ll| and later 
extended by Kluyver (1906) [l2| and Rayleigh (1919) 
it is assumed that particles perform random displace- 
ments ('flights') along straight lines, and that at the end 
of each flight (a 'collision' with the surrounding medium) 
the direction of propagation changes at random. 

When the number of transported particles is much 
smaller than the number of the particles of the interacting 
medium, so that inter-particles collisions can be safely ne- 
glected, it is reasonable to assume that the probability of 
interacting with the medium is Poissonian. For the case 
of neutrons in a nuclear reactor, e.g., the ratio between 
the number of transported particles and the number of 
interacting nuclei in a typical fuel/moderator configura- 
tion is of the order of 10~^^, even for high fiux reac- 
tors 3 . It follows that flight lengths between subsequent 
collisions are exponentially distributed (hence we will call 
this process exponential flights in the following [l^l)- We 



assume that collisions can be either of scattering or ab- 
sorption type. At each scattering collision, the flight di- 
rection changes at random, whereas at absorption events 
the particle disappears and the flight terminates. Each 
flight can be seen as a random walk in the phase space 
of position r and direction w in a d-dimensional setup. 

The particle density 5'(r,a;,<) represents the probabil- 
ity density of finding a transported particle at position 
r having direction w at a time t, up to an appropriate 
normalization factor. In many applications, the actual 
physical observable is the average of the density ^'(r, w, t) 
over the directions u, namely 



*(r,t) 



1 



'^{v,uj,t)duj. 
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where f7<j = / = 27r''/Vr (d/2) is a normalization fac- 
tor corresponding to the surface of the unit d-dimensional 
sphere and r() is the Gamma function [l4| . 

Along with the development of Monte Carlo methods, 
numerical solutions 4'(r,<) to complex three-dimensional 
linear/nonlinear transport problems coming from applied 
sciences are becoming accessible to an high de gree of ac- 
curacy: criticality calculations in reactor cores [la|, scat- 
tering and absorption in heated plasmas [165, propaga- 
tion through anisotropic scattering centers in atmosphere 
or fluids [l7[ ■, and charge transport in semiconductors un- 
der external fields |18l - [2(]| . only to name a few. Nonethe- 
less, even for the simplest systems, many theoretical 
questions remain without an answer, so that the study of 
exponential flights has attracted a renovated interest in 
recent years; see, e.g., f2T| - [23 |. In particular, it has been 
emphasized that the dimension d deeply affects the na- 
ture of ^{t, t), and prevents in most cases from obtaining 
explicit results. The aim of our work is to investigate ex- 
ponential flights in a generic d-dimensional setup, under 
simplifying hypotheses. Here, we will mostly focus on 
establishing insightful relationships between space, time 
and the statistics of particle collisions within a given vol- 
ume. A number of new results will be derived, concerning 
unbounded, bounded, scattering as well as absorbing do- 
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FIG. 1: The free propagator 5'(2:|n) in the transformed space. 
Left: for increasing values of the dimension, d = 1 

(blue), 2 (red), 3 (green), and 4 (black). Right: ■^{z\n) for 
d = 3 and increasing number of collisions, n — 1 (blue), 2 
(red), 3 (green), and 10 (black). 



This paper is structured as follows. In Sec. |TT] we re- 
call the mathematical formalism, introduce the physical 
variables and derive their inter-dependence for any d. In 
Sec, mil we detail the structure of the spatial moments of 
the particle ensemble. Sec. IIVI is devoted to the analy- 
sis of the collisions statistics in a given domain. Then, 
in Sec. |V]we examine the distinct cases d = 1,2,3 and 
4. Both the spatial and temporal evolution of the par- 
ticle ensemble are considered, and results for bounded 
domains are obtained by resorting to the method of im- 
ages. We provide a comparison between analytical (or 
asymptotic) findings and Monte Carlo simulations. Con- 
clusions will be finally drawn in Sec. IVII 



II. GENERAL SETUP 

Within the natural framework of statistical mechanics, 
the evolution of the particle density 5'(r,a;,t) for expo- 
nential flights is governed by the linear Boltzmann equa- 
tion 0. Linearity stems from neglecting inter-particle 
collisions. In the hypothesis that an average particle en- 
ergy can be defined (the so called one-speed transport), 
and that the physical properties of the medium do not 
depend on position nor timejthe Boltzmann equation for 
the density 5'(r,w,f) reads %\2)\ 

1 d 

- — ^'(r,tj,t)+a;- Vr*(r,a;,t) = 
V at 

= -^(^(r,^,^) -l-cr j k{uj' ,bj)-^{r,bj\t)dbj' + (2) 

where at is total cross section of the traversed medium 
(carrying the units of an inverse length), a is the scat- 
tering cross section, v is the particle speed, and S is the 
source. The total cross section <Tt is such that l/ut rep- 
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FIG. 2: The free propagator ^(r\n) for d = 1, with n = 
1 (blue), 5 (red), 10 (green), and 50 (black). Monte Carlo 
simulation results are displayed as symbols. The solid line 
is the theoretical result, Eq. ((50]). The dashed line is the 
diffusion limit, Eq. (fTU)) . 



resents the average flight length between consecutive col- 
lisions (the so-called mean free path), and is related to 
the scattering cross section a and to the absorption cross 
section (Ja hy at = a + aa- The quantity k{ui' ,ui) is the 
scattering kernel, i.e., the probability density that at each 
scattering collision the random direction changes from w' 

to UJ. 

Denoting by ^'(r,w,t) the solution of the Boltzmann 
equation ^ for a medium without absorptions [a a — 0), 
the complete solution with absorption ^'^(r,^;,^) can be 
easily obtained by letting 



*a(r,t^,i) = «'(r,tj,t)e" 



(3) 



thanks to linearity [21]. This allows primarily addressing 
a purely scattering medium {at — a), without loss of 
generality. 

At long times, i.e., far from the source, the direction- 
averaged particle density ^'(r, t) is known to converge to 
a Gaussian shape, namely. 



*(r,t) 



(47rDt)" 



(4) 



the quantity D = v/ {da) playing the role of a diffusion 
coefficient 8]. However, Eq. (j4]) is approximately valid 
for rcr 3> 1, and can not capture the particle evolution at 
early times, nor the finite-speed propagation effects. In- 
deed, diffusion implicitly assumes a non- vanishing prob- 
ability of finding the particles at arbitrary distance from 
the source. Deviations from the limit Gaussian behav- 
ior are well known, e.g., for neutron Q as well electron 
transport [l9| . 

In the following, we outline the relation between 4'(r, t) 
and the underlying exponential fiight process. 
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A. The free propagator without absorptions 

Consider a c?-dimensional setup. A particle, originally 
located at position Tq in a given domain, travels along 
straight lines at constant speed v, until it collides with the 
medium. The position of a particle at the n-th collision 
can be expressed as a random walk r„ = Tq + 
i.e., a sum of random variables r^. The flight lengths 
£ = \ri — ri_i| are assumedly identically distributed and 
obey an exponential probability density 



ate 



(5) 



with CTt > 0. The exponential law in Eq. ([5]) stems 
from assuming a uniform distribution of the scattering 
centers in the traversed medium. Heterogeneous materi- 
als, such as complex fluids, would generally lead to clus- 
tered scattering centers, obeying to, e.g., negative bino- 
mial distributions, and in turn to non-exponential flight 
lengths However, we will focus our attention on 

homogeneous media. 

At each collision, the particle randomly changes its di- 
rection according to the scattering kernel k{uj',Lu). For 
the sake of simplicity, we assume here that the scattering 
is isotropic, so that k{uj',uj) has a uniform distribution, 
independent of the incident direction u' . 

Once a flight length has been sampled from (p{£), the 
new direction uj is therefore uniformly distributed on 
the d-sphere of surface i'^~^ftd- Therefore, by virtue of 
the apparent spherical symmetry, the transition kernel, 
i.e., the probability density of performing a displacement 
from ri_i to r^, depends only oni — jr; — ri_i|, and reads 



(6) 



We initially neglect absorptions, so that at = <t: in 
one-speed transport, this condition can be either seen 
as the particles been scattered, or equivalently being ab- 
sorbed and then re-emitted (with the same speed) at each 
collision. This latter interpretation would correspond, 
e.g., to a criticality condition in multiplicative systems 
for neutron transport. 

We define then the free propagator ^'(r|n) as the prob- 
ability density of finding a particle at position r at the 
ri-th collision, for an infinite medium, i.e., in absence of 
boundaries. We adopt here the convention that the parti- 
cle position and direction refer to the physical properties 
before entering the collision; for instance, the index n = 1 
refers to uncollided particles, i.e., particles coming from 
the source and entering their first collision. 

Assuming that all the particles are isotropically emit- 
ted at ro — 0, the particle density 4'(r|n) must depend 
only on the variable r — |r|, because of the spherical sym- 
metry. On the basis of the properties exposed above, the 
particle propagation as a function of the number of col- 
lisions is a Markovian process in the variable r„, where 
for each collision i = 1, n the new propagator is given 




FIG. 3: The free propagator '^(r\n) for d = 2, with n = 
1 (blue), 5 (red), 10 (green), and 50 (black). Monte Carlo 
simulation results are displayed in symbols. The solid line 
is the theoretical result, Eq. (|62p . The dashed line is the 
diffusion limit, Eq. (fTU)) . 



by the convolution integral 

*(r|i) = J TTrfdr - r'|)*(r'|i ~ l)dr' , (7) 

with initial condition ^(r|0) — 6{r). In particular, by di- 
rect integration we immediately get the uncollided prop- 
agator 



«'(r|l) = TTd{r). 



(8) 



It is convenient to introduce the d-dimensional Fourier 
transform of spherical-symmetrical functions, as in the 
subsequent analysis this will allow easier deriving of 
the properties of the exponential flights. Denoting by 
z the transformed variable with respect to r, for any 
spherical-symmetrical function /(r) we have the follow- 
ing transform and anti-transform pair f(z) = {/('")} 
and fir) = {/(z)} i 

+ 00 



/(z) = z^-^l^ {2nf^ / T'"\h/2-i{zr)f{r)dT 
Jo 

fir) = (27r)"'^/2 / z'^/V,/2_i(rz)/(z)dz, (9) 



where Ji,() is the modified Bessel function of the first 
kind, with index v [14| . It is apparent from Eqs. (|9]) 
that the dimension d of the system plays a fundamental 
role, in that it affects both the transition kernel and the 
Fourier transform kernel itself. 

The convolution integral in Eq. ([7]) in Fourier space 
gives the algebraic relation 



^iz\i) = 7rd(z)*(z|i- 1), 



(10) 
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i > 1, with initial condition 4'(z|0) = 1. By recursion, it 
follows that in the transformed space 



*(z|n) = TTdizY 



(11) 



It turns out that the Fourier transform of the transition 
kernel TTd{z) can be explicitly performed in arbitrary di- 
mension, and gives 



2' 



(12) 



where 2^1 is the Gauss hypergeometric function [T3 |. 
Hence, we finally obtain 



^-(21^1) 



1 d z^ 
"i^i 7:, 1, 7:; 7 



(13) 
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The quantity 'Kdiz) is positive for d > 1; moreover, 
■nd{z = 0) = 1, which ensures normalization and posi- 
tivity of the propagator. In Fig. [T] we visually represent 
the effects of dimension and number of collisions on the 
shape of ^'(z|?i). Remark in particular that the spread 
of 5'(z|n) increases with d, for a given n. On the con- 
trary, 4'(z|n) becomes more peaked around the origin 
with growing n, for a given d. 

Formally, performing the inverse Fourier transform of 
Eq. ([T^ gives the propagator '^{r\n) = {'^{z\n)} 
for an arbitrary d-dimensional setup. Actually, in some 
cases this task turns out to be non-trivial. Nonetheless, 
even in absence of an explicit functional form for the 
propagator, information can be extracted by resorting 
to the Tauberian theorems. In particular, the expansion 
of '^{z\n) for z/a ^ 1 gives the behavior of ^'(r|n) for 
rcr ^ 1, i.e., far from the source, in the so-called diffusion 
limit 0, viceversa, the expansion of ^(z|n) for z/cr ^ 
1 gives the behavior of 5'(r|ri) for ra <^ 1, i.e., close to 
the source. We recall that 2F1 is defined through the 
series [iJI 



At the leading order for z/a ^ we therefore have 



(14) 



(15) 



Observe that Eq. ([T5| can be viewed as the expansion 
of an exponential function. Then, the inverse Fourier 
transform gives the Gaussian shape 



FIG. 4: The free propagator '^(r\n) for d = 3, with n = 
1 (bhie), 5 (red), 10 (green), and 50 (black). Monte Carlo 
simulation results are displayed in symbols. The dashed line is 
the diffusion limit, Eq. (I74p . The solid line is the approximate 
propagator, Eq. (f7D)) . 



attraction of the Central Limit Theorem Q . We mention 
that clustered scattering centers, with non-exponential 
flight lengths, may lead to non-Gaussian limiting statis- 
tics [17[. Remark the close analogy between Eqs. 



and (|16p : in particular, V and D differ by a factor av, 
which roughly speaking represents the average number of 
collisions per unit time. 

Moreover, at the leading order for z/cr —> 00 we have 
the expansion 



7rrf(z) 



0^r(f) 
r(^) 



(5). (2-.) (5)%..., (17) 



where the first term vanishes for d ~ I. 
Fourier transforming, we have for ru <C 1 



^{r\n) 



71. d 



[ray 



By inverse 



(18) 



when > 1, and ^{r\n) 



n.l / \ 

C2 [rcr) 



2n-l 



when 

d = 1. Here c"''^ and are constants depending on n 
and d. It can be shown that the divergence at the origin 
in ^(r|n) due to the Dirac delta source disappears after 
n > d collisions for d> 1, and n > 1 for d = 1. 



Relation between collision number and time 



^{r\n) 



(AnVn) 



d/2 ' 



(16) 



which is valid for rcr ^ 1,1) — l/{da'^) playing the role of 
a diffusion coefficient. This stems from the exponential 
flights having finite-variance increments, < -t-00, so 
that their probability density 4'(r|n) falls in the basin of 



Assume again that at = cr, i.e., that there are no ab- 
sorptions. The free propagator ^(r|n) gives information 
on the position of a transported particle at the moment of 
entering the n-th scattering collision. The link between 
the travelled distance, the flight time and the number of 
collisions is provided by the speed v. Indeed, once a flight 
of length £ between any two collisions has been sampled 
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FIG. 5: The collision density ^'(r) for d = 3, with increasing 
number of summed collisions, A'' = 10^ (blue), 10"^ (red), 10'* 
(green), and 10'' (black). The dashed lines are the asymptotic 
limits close (Eq. ([79])) and far (Eq. ([78])) from the source. 



from ip{£), the flight time must satisfy ti = ti—ti-i = Ijv. 
Hence, the transition kernel, i.e., the probability density 
of performing a displacement from r^^i at ti^\ to at 
ii, will be given by 



(19) 



It follows that inter-collision times are exponentially dis- 
tributed 



(20) 



where r — \/{(jv) represents the average time between 
collisions. 

We define the propagator ^'(r, t|n) as the probability 
density of finding a particle at position r at time t at the 
n-th collision. From the Markov property of the process, 
at each collision i — 1 , . . . , n we have 

*(r,t|i) = j dT' j dt'Tid{\v - T'lt - t')-^{r\t'\i ^ I), 

(21) 

with initial condition *(r, t|0) = 5{r)5{t/T) = TS{r,t). 
In particular, by direct integration we immediately get 
the uncollided propagator 

^'(r, t\l) = 7rd(r, t) = T*(r|l),5 (i - ^ • (22) 

We denote the Laplace transform of a function g{t) by 
its argument, i.e.. 



r + co 

g[s)^C{g{t)}^ e-^'g 
Jo 



{t)dt. (23) 



Then, from the double convolution integral in Eq. ((2T|) 
we have the algebraic product in the Fourier and Laplace 
space 



*(z,s|i) = 7rd(z,s)*(z,s|i- 1), 



(24) 



i > 1, with ^'(z, s|0) = r. By recursion, it follows that in 
the transformed space we have 



^'(2,s|n) = Tiid{z,sY' 



(25) 



It turns out that the Fourier and Laplace transform of 
the transition kernel 7r(j(z, s) can be explicitly performed 
in arbitrary dimension, and gives 



2Fi(i,l,'^- 



l + ST 

where Q{s) = a{\ + sr). Hence, we finally obtain 



(26) 



^'(z,s|n) 



2F1 



(1 1 d.. 

' ' 2 ' 



\ + ST 



(27) 



Moreover, the following relation follows: ^'(2,5 
0|n) = T^{z\n), so that 



'^{r\n) 



+00 



^{r,t\n)dt, 



(28) 



i.e., '^{r\n) can be intepreted as the time average of 
^(r, t|n). Finally, the propagator 5'(r, t) will be given by 
the sum of the contributions 4'(r, t|n) at each collision, 
namely 



«'(r,<) = ^*(r,t|n). 



(29) 



n = l 



Taking the Fourier and Laplace transform of Eq. (|29l) . 
we then have 



*(z,s) = ^^'(z,s|n) = T- 



i^d{z,s) 



(30) 



with *(r, t) = C-^T^' {^{z, s)} 



C. Absorptions 

In presence of absorptions {a a > 0), the propagator 
^'a (f , t) can be obtained from Eq. ([3]) by integrating over 
directions. This relation holds true at each collision, so 
that we have 



v|/^(r, t|n) = ^{r,t\n)e 



-tjr 



(31) 



where Ta = l/{(Tav). Hence, from Eq. ((28l) . by replacing 
a with at we get for the propagator 5'a(z|n) 

1 /■+°° 1 1 

^aiz\n) = — / ^'aiz,t\n)dt= —^'{z,s= —\n), 

n Jo Tt Ta 

(32) 
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where tj = 1/ {utv) represents the average flight time be- 
tween any two colhsions. Now, observe that from Eq. (P7)) 
we have 



* ( z,s = —\n 

Ta 



T* I z — \n 



(33) 



where the quantity p — a/at, < p < I, can be inter- 
preted as the the probabihty of being scattered, i.e., not 
being absorbed, at any given cohision. Then, by remark- 
ing that t/ti — 1/p, it follows that 



*a(-z|n) = ^'(pz|n)p" 



(34) 



The propagator in Eq. (|34l) is then given by the prod- 
uct of the free propagator, with the total cross section at 
replacing the scattering cross section cr, times the prob- 
ability of having survived up to entering the n-th colli- 
sion. Remark that the total cross section and the non- 
absorption probability are related by cr = pat- When 
the absorption length is infinite, (Xa — 0, p — 1 and 
we recover the free propagator for pure scattering, with 
at = a. 



D. Boundary conditions 

So far, we have assumed that the medium where par- 
ticles are transported has an infinite extension, hence 
the name free propagator. More realistically, we might 
consider finite-extension media with volume V enclosing 
the source, so that boundary conditions come into play 
and affect the nature of the propagator. Several bound- 
ary conditions can be conceived according to the spe- 
cific physical system, among which the most common 
are reflecting and leakage. This issue has been exten- 
sively examined, e.g., for radiation shielding in Reac- 
tor Physics and for electron motion in semiconduc- 
tors [2^, HBl- Here we will focus on leakage boundary 
conditions, where particles are considered lost as soon as 
their trajectory has traversed the outer boundary dV of 
the domain. While the volume V is in principle totally 
arbitrary, in the subsequent calculations for the sake of 
convenience we will assume that V = V{R) is a sphere of 
radius R centered in Tq = 0. 

From the point of view of the propagator, leakages 
can be taken into account by assuming that the popu- 
lation density ^(r|n) at any n vanishes at the so-called 
extrapolation length rg, i.e., ^{r = re|n) = 
Because trajectory do not terminate at the boundary, 
but rather at the first collision occurred outside the vol- 
ume, the extrapolation length is expected to be larger 
than the physical boundary of V{R) and can be deter- 
mined from solving the so-called Milne's problem associ- 
ated to the volume 27, 28]. In general, re is of the kind 

= R[\ + Udl [R'y)], the dimensionless constant > 
depending on the dimension of the system [28| . When the 
scattering length is much smaller than the typical size of 
the volume, i.e., aR ^ 1, the extrapolation length coin- 
cides with the physical boundary, Ve — > R- This means 



that the inter-collision length is so small as compared to 
the total travelled distance that the first collision out- 
side the domain will actually be very close to the last 
collision inside the domain, which corresponds to ^(r|n) 
vanishing at R 



III. SPATIAL MOMENTS OF THE 
PROPAGATOR 

The moments of the propagator provide an estimate 
of the spatial evolution of the particles ensemble, as a 
function of the number of collisions or time. Due to the 
supposed spherical symmetry, we expect all the odd mo- 
ments to vanish. We define the m-th moment of fir) 
over the spherical shell r'^^^H.ddr as 



r 



m-\-d- 



'f{r)dr. 



(35) 



Performing the integral in Eq. p5p would require explic- 
itly knowing f{r). However, as shown in the Appendix [XI 
the TO-th spatial moment /(r) of can be expressed as a 
function of the m-th derivative of f{z) with respect to z 



(r™) = 



0Fr(^) 



r(f)r(i±2i) d(iz) 



(36) 



when m is even, and (r™) — otherwise. By setting 
f{z) = '^{pz\n)p"~^ in Eq. ([36|) we then have {r"^){n) as 
a function of the number of collisions. Analogously, by 
setting f{z) = ■^aiz,s) = ^{z,s + l/r^) we get (r'")(s), 
which gives the evolution as a function of time, upon in- 
verse Laplace transforming. In particular, for the spread 
TO = 2 of the propagator with absorptions we get 



{r^)in) = -,p- 



and 



,-t/r _ 1 



(37) 



(38) 



In absence of absorption, Ta 00 and p — >■ 1, the par- 
ticle spread {r'^){n) is linear with respect to n. On the 
contrary, {r'^){t) has a ballistic behavior (c>c t'^) at earlier 
times (where transport is dominated by velocity), and a 
diffusive behavior (cx t) at later times (where transport 
is dominated by scattering). The transition between the 
two regimes is imposed by the time scale r. A remark- 
able feature is that in either case the spread does not 
explicitly depend on the dimension d. Intuitively, this 
can be understood by considering that (independently of 
the dimension d of the embedding setup) the vectors 
and ri_|_i define a plane with random orientation, so that 
space is explored by plane surfaces at each collision. This 
is in analogy, e.g., with the behavior of d-dimensional 
Brownian Motion 
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IV. COLLISION DENSITY AND COLLISION 
STATISTICS 

In many physical problems, one is interested in assess- 
ing the statistics of the time ty or the collision number ny 
spent inside a given domain V. A prominent role in char- 
acterizing a physical system is played in particular by the 
mean residence time (iy) and the mean collision number 
(riy) [29]. In Reactor Physics, for instance, the average 
number of neutron collisions within a region would be 
related to such issues as the nuclear heating or nuclear 
damage in fissile as well as structural materials [1, 0] ■ We 
introduce the collision density ^'(r) JscI], which is defined 
as 



to the formulation of first-passage problems, i.e., deter- 
mining the distribution of the time, or collision number, 
required to first reach the boundary. Consequently, (ty) 
is preferentially called mean first-passage time . 

It can be shown that higher order moments of ny can 
be obtained by recursion, in terms of Kac integrals analo- 
gous to those in Eq. (|4T|) . A detailed treatment is beyond 
the scope of the present paper and is discussed in [33] ■ 

The previous discussion shows that 4'(r) is key in de- 
termining residence times and collisions statistics. After 
formally carrying out the summation in Eq. ([39]) , we have 
= Trd{z)/[1 — TTd{z)], so that \E'(r) for a free propa- 
gator in absence of absorption is defined in terms of the 
following Fourier integral 



N 

vE'(r) = lim "iiMn). 

n=l 



(39) 



From Eq. (p9)) . it follows that ^'(r) can be equivalently 
obtained from the propagator 5'(r, i) as 



lim — 



*(r,i)di. 



(40) 



For an infinite 'observation time' T, the mean residence 
time (tv) within V can be expressed in terms of the colli- 
sion density \l/(r) in the same domain by slightly adapting 
an argument due to Kac [3l|, [12] 



(41) 



where ri = ]ri — roj. 

As for the collision number, the probability of perform- 
ing ny collisions in V is given by 



Vinv)^ / dr*(rlny)- / dr^{r\nv + l), (42) 
Jv Jv 



hence the moments 



(nv) 



+ 00 

E 



^Vinv) 



(43) 



(2^) 



d/2 



( \ 

z'"'j,/2-i{rz) , ""'^"l , dz, (45) 

Trd[z) 



1 



whose convergence depends on the dimension d of the sys- 
tem. It turns out that convergence is ensured for d > 2, 
which means that for Id and 2d finite-size domains with 
transparent boundaries {nv) and {tv) diverge as — > oo. 
This result is in analogy with Polya's theorem, which 
states that random walks on Euclidean lattices are re- 
current for d < 2 0. As shown in the following, we can 
nonetheless provide an estimate of such divergence as a 
function of A'^, i.e., single out a singular term from a func- 
tional form. For finite domains with leakage boundary 
conditions and/or absorptions {a a > 0), ^'(r) is defined 
also for Id and 2d systems. For d > 2, Tauberian the- 
orems show that the asymptotic behavior of Eq. (j^S]) is 
given by 



for large r, and 



*(r) 



£(|) 



{ra 



,2-d 



(46) 



(47) 



close to the source. 



From the definition of 5'(r), it follows immediately that 



{nv) 



{tv) 



(44) 



i.e., the integral of the collision density over a volume V 
gives the mean number of collisions within that domain, 
hence the name given to ^'(r). 

Remark that both (ty) and {nv) depend on the bound- 
ary conditions imposed on dV, which in turn affect the 
functional form of the propagator, and then ^(r). Using 
a free propagator corresponds to defining a fictitious vol- 
ume V, whose boundaries dV are 'transparent', so that 
particles can indefinitely cross dV back an forth. On the 
contrary, the use of leakage boundary conditions leads 



V. ANALYSIS OF d-DIMENSIONAL SETUPS 

In the following, we detail the results pertaining to 
specific values of d. We choose l/at as length scale and 
we work with dimensionless spatial variables r — rut. 
Remark that in absence of absorption the length scale is 
1/ct, since p = 1- 



A. One-dimensional setup, d = 1 

The case d = 1 allows illustrating the general struc- 
ture of the calculations. One potential application of 
this framework could be provided by nanowires or carbon 
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nanotubes (almost Id systems) in electron transport p^ . 
The transition kernel is 



whose Fourier transform is 



7ri(z) = 



l + z2 



(48) 



(49) 



Starting from ^(z, 0) = 1, the free propagator ^(r|n) can 
be explicitly obtained by performing the inverse Fourier 
transform of 'i'{z\n) — 7ri(z)", and reads 



22- 



v^r(n) 



(50) 



where Ki,{) is the modified Bessel function of the sec- 
ond kind, with index 



The same formula has 
been recently derived, e.g., in [^J] as a particular case 
of a broader class of random flights. In Fig. [2] we pro- 
vide a comparison between Monte Carlo simulation re- 
sults (symbols), the analytical formula Eq. (150)) (solid 
lines), and the diffusion limit, Eq. (fT6|) (dashed lines), 
for different values of n. Remark in particular that the 
diffusion limit is not accurate for small n, and becomes 
progressively closer to the exact result for increasing n, 
as expected. At intermediate n, the tails of the propa- 
gator ([50)1 are always fatter than those predicted by the 
diffusion approximation. 

In a Id setup, the collision density '${r) for the free 
propagator diverges. Nonetheless, it is possible to single 
out the divergence as follows 



N 

\l>(z) = lim \l>(z|n) 

n=l 



-N 



lim 5^ 5-^^ . (51) 



For fixed N, the inverse transform can be explicitly per- 
formed in terms of hypergeometric functions. Retaining 
the non- vanishing terms for large N, we have 



r(i + ^) r 
r{N)^ 2 



N r 
¥ " 2' 



(52) 



which is composed of a term diverging with V-^V (not 
depending on r), and a functional part which is linear in 
r (not depending on N). 

For the propagator with absorptions, from Eq. (I34p we 
have 



AT 



■^a{z) = lim y"-^{z\n)f'-^ 



n=l 



1 — p + 



(53) 



Then, performing the inverse Fourier transform, we get 



*a(r) 



/I— pr 



(54) 



Remark that Eq. ((M)) has been derived, e.g., in [34| by 
solving the stationary Boltzmann equation in Id. When 



p — ^ 0, the particles are almost surely absorbed at the 
first collision, and we have the expansion 



T / \ e /-, P pi" 



(55) 



so that at the first order the collision density has the 
same functional form as the uncollided propagator. At 
the opposite, when p — )■ 1 the particles are almost surely 
always scattered (ct — ^ cr), and we have the expansion 

^air)^—^-l + ..., (56) 
2vl-p 2 

and '^a{r) diverges as the collision density associated to 
the free propagator, as expected. 

The case of leakage boundary conditions can be dealt 
with by imposing that the propagator 4'(r|n) must vanish 
for any n at the extrapolated boundary r^. For d — 1, 
the extrapolated length is given by re — R[l + Ud/R] 
with ui — 1 (28| . i.e., the propagator must vanish at one 
scattering length outside the physical border R of the do- 
main. By resorting to the method of images Q, which 
allows solving for the propagator in presence of bound- 
aries in terms of the propagator in absence of boundaries, 
we therefore have for the collision density 



*fl(0 



(57) 



for r < R, and 5'(r) =0 elsewhere. 

The moments of the residence time (or, equivalently, of 
the number of collisions) within a sphere of radius R can 
be explicitly computed based on Eq. (j44p . The moments 
associated to the free propagator clearly diverge. Here 
we separately analyze the propagator with absorptions 
(in an infinite domain) and the propagator with leakages 
at the boundary r — > R (without absorption). For 
the case of absorptions, the average residence time within 
a (fictitious) sphere of radius R reads 



l-p 



-Tt, 



(58) 



assuming that particles can cross the boundaries of the 
sphere an infinite number of times. When the radius of 
the sphere is large as compared to the typical particle 
displacement, we have {tji) ~ 1/(1 — p)Tt, which gives 
(tfl) ~ Tt when p — ^ 0, and diverges for p — > 1. For the 
case of leakages, the mean first-passage time reads 



R{2re-R) 



(59) 



When the radius of the sphere is large as compared to 
the typical particle displacement, — >■ R, and we have 



B. Two-dimensional setup, d = 2 

The case d = 2 has a key interest in assessing, e.g., 
the dynamics of chemical and biological species on sur- 
faces [9|. Moreover, it concerns also quantum wells and 
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inversion layers in electron transport fioj . The transition 
kernel reads 



whose Fourier transform is 



2ttV 



(60) 



where Te = + U2/-R] and the Milne's constant U2 — 
l-2/7r2 m. 

The moments associated to the free propagator clearly 
diverge. For the propagator with leakages at the bound- 
ary r = Tf, > R the mean first-passage time within a 
sphere of radius R reads 



(61) 



Starting from ^'(z, 0) = 1, the free propagator '^(r\n) can 
be explicitly obtained by performing the inverse Fourier 
transform of '^{z\n) = 7r2(z)", and reads 



^^(rln) = 



2-f r-i+t (r) 
TrFf^i) 



(62) 



This result was previously established by [35[, and has 
later appeared in, e.g., [H, H^. In Fig. [3] we compare 
the Monte Carlo simulation results (symbols) with the 
theoretical formula in Eq. ([62]) . for different values of n. 
The diffusion limit, Eq. (jl6p . is also shown in dashed 
lines. Remark that the diffusion limit is not accurate for 
small 71, and becomes progressively closer to the exact 
result for increasing n. At intermediate n, the tails of the 
propagator (p^ are always fatter than those predicted by 
the diffusion approximation. 

In a 2d setup, the collision density \E'(r) diverges. 
Nonetheless, analogously as done for the Id case, it is 
possible to single out the divergence as follows 



N 



^'(z) = lim ^'(z|n) = lim 



1-1 



-N/2 



- 1 



(63) 

For fixed N, the inverse transform can be explicitly per- 
formed. Details of the rather cumbersome calculations 
are left to the Appendix [BJ Retaining the non- vanishing 
terms for large N, we have 



*(r) 



log(iV) e-'' Ei (-r) - 21og(r) 
27r 27rr 27r 



(64) 



where Ei is the exponential integral function To 
the authors' best knowledge, the formula for the colli- 
sion density \I'(r) has not appeared before, and might 
then provide a useful tool for describing the migration of 
species on 2d environments. Similarly as in the Id case, 
^(r) is composed of a term diverging with log N (not de- 
pending on r) , and a functional part in r (not depending 
on TV). In deriving Eq. (|M)) we have neglected a con- 
stant term which is small compared to log(A^), namely, 
[log(2) — 7/2]/7r, where 7 ~ 0.57721 is the Euler's gamma 
constant [14| . 

The collision density with leakages at r — can be 
obtained again by the method of images, whence 



*i?(0 



e'' 
2'Kr 



Ei (-r) -Ei(- 



21og(^) 



2Txr„ 



27r 



(65) 



1 - e-" + Re-^^ + R^ - 



{tn) = 



R2E^ (-R) - R^Ez {-r,) - 2R' log(ii) 



r. (66) 



When i? >• 1, we have — R, and we get {tfi} ~ (1 -I- 
R^)t/2 in the diffusion limit. 

As for the collision density with absorptions, calcula- 
tions analogous to the Id case lead to 



*a(r) = ^ ^ + / dz. 

TT 27r 7o VTT^ + p 

(67) 

We could not find an explicit expression for the latter 
integral in terms of elementary functions. However, the 
limits for small and large scattering probability can be 
easily obtained, and read 



2nr TT 



(68) 



when p — > 0, and 



27r 27rr 



2tt 



when p — >■ 1, respectively. The former expression gives 
the uncollided propagator at the leading order, whereas 
the latter diverges logarithmically as p ^ 1. 



C. Three-dimensional setup, d = 3 

The case d = 3 plays a prominent role, among others, 
in Reactor Physics, in that it concerns the transport of 
neutrons and photons through matter [6;-^] , and is key in 
describing electron transport in bulk semiconductors [13] ■ 
On the basis of the strikingly similar form of Eqs. ([50)1 
and (1621) , it would be tempting to postulate an analogous 
expression for the propagator in 3d. For d = 1 we have 
indeed the functional form 4'(r|n) cx r~^/^~^"-K_i/2+n(.i^)^ 
and for d = 2 *(r|n) cx r-^+"/'^K_i+n/2{r). Then we 
could conjecture an exponent — d/2 + n/d, so that 



v['*(r|n) 



-^/^+"j^-i/2+n(r-) 
25+t7rir(?) 



(70) 



by imposing normalization. Unfortunately, this is not the 
case, and ^'*(r|n) is not the true 3d propagator. Actu- 
ally, few explicit results can be derived, and much of the 
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analysis is therefore devoted to the asymptotic behavior. 
The transition kernel reads 



whose Fourier transform is 
7^3(2) 



arctan (z) 



(71) 



(72) 



The propagator 5'(r|n) with initial condition ^(z, 0) = 1 
involves then the following integral 



^'(r|n) 



1 



27r2r 



+00 



z sin (rz) 



arctan (z) 



dz, (73) 



which can not be carried out explicitly for arbitrary n. 
In the diffusion limit z ^ 1 we have [arctan(z)/z]"' ~ 
1 — nz^ so that 



*(r|n) 



8n3/27r3/2 ' 



(74) 



as expected from Eq. pB)) . In Fig. S] we compare the 
Monte Carlo simulation results (symbols) with the diffu- 
sion limit, Eq. (|74|) (dashed lines), and with the approx- 
imate propagator, Eq. ([70)) (solid lines). Remark that 
Eq. ([70]) provides a fairly accurate approximation of the 
simulation results, except close to the source. 

After carrying out the sum over n, the collision density 
^'(r) is given by the following integral 



l'(r) 



1 



. , , arctan (z) 
, zsm(rz) ■ r^az. 







z — arctan (z) 



(75) 



which again can not be performed explicitly [30| • As be- 
fore, we consider then the asymptotic behavior. Denoting 
h{z) = arctan(z)/[z — arctan(z)], we have 



hiz) 



36 
175 



in the diffusion limit z <C 1, and 



h{z) 



2^ 



4z2 



8zJ 



(76) 



(77) 



close to the source. Similar expansions appear, e.g., 
in [soj . as derived from the analysis of the Boltzmann 
equation. Correspondingly, by performing the respective 
integrations we have 



*(r) ~ — 
for r 3> 1, i.e., far from the source, and 



(78) 



for r ^ 1, i.e., close to the source. Remark that for 
d = 3 4'(r) does not diverge even for infinite domains 
without absorptions. In Fig. [S] we compare the Monte 
Carlo simulation results (symbols) with the asymptotic 
limits close to and far from the source, Eqs. (|79)) and ([78]), 
respectively (dashed lines). The simulation results pro- 
gressively approach the asymptotic limits as the number 
N of summed collisions increases. 

Eqs. ([7^ and ([75]) can provide asymptotic estimates 
for the collision density with leakages at the boundary 
r = R. By the method of images, we have that the col- 
lision density with boundaries is ^_R(r) = \['(r) — 5'(re), 
with re = R[l + U3/R] and U3 - 0.7104 [13 . 

The moments of the residence time within a sphere of 
radius R can be explicitly computed based on Eq. ([^1)1 
for the free propagator \I'(r|n), i.e., when particles can 
freely cross the surface of the sphere. We have 



when i? ^ 1, and 



R 



-R']t 



(80) 



(81) 



for i? <C 1. Moreover, for leakage boundary conditions 
at the surface, from 5'/j(r|n) we have 



i?2 /3 



when i? ^ 1, and 



for i? < 1. 



,7r2-4 i?2re(4-7r2)-4 



8 
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(82) 



T (83) 



D. Four-dimensional setup, d = 4 

The case d = 4 is briefly presented here for the sake of 
completeness. The transition kernel reads 



whose Fourier transform is 
7r4(z) 



27r2£3 ' 



l + VT 



(84) 



(85) 



We could not find an explicit representation for the in- 
verse Fourier transform of ^'(z|n) — TT4{z)". Nonetheless, 
the propagator ^'(r, tin) is known and reads 
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for vt > r [2i|. Hence, it follows that the propagator 
^{r\n) = J ^{r,t\n)dt/T can be obtamed from solvmg 
the integral 



'i'{r\n) 



7r2r {n - 1) 



+ 00 



e-'z-'-'Uz^-r^) dz. 



This integral can be performed, and gives 

1 



(87) 



*(r|n) 



247r3/2 



[A + B + C], 



where 



B 



4„n-4 



-2V 



C = 22 



r(2-n) 



n ]_ -2+7t . 
' 2 ' 2 ' 2 '4 



1^5 



/ 1-n 
V 2 ' 



— n 3 
2 ' 



-1+n . r" 



(89) 



1F2O being an hypergeometric function 14| . 
As for the collision density ^'(r), we have 



1 



Ji (rz) (1 + Vl + z^) (90) 



27r2r JO 

which can be computed explicitly and gives 

e""" 1 

^(r) = 



27r2r3 



(91) 



Finally, the collision density in presence of leakages at 
the boundary r = R can be obtained by resorting to the 
method of images, whence \E'_R(r) = '^{r) — '^{re), with 
Te = R[l + U4/ R]. The constant U4 has been estimated 
by running a Monte Carlo simulation and determining 
the extrapolation length, and reads U4 ~ 0.5. 

The moments of the residence time within a sphere of 
radius R can be explicitly computed based on Eq. pi)) 
for the free propagator ^(r|n), namely. 



(tn) = {1 + R'- e-^) 



(92) 



Moreover, for leakage boundary conditions at the surface, 
from ^ij(r|n) we have 



l + R^ 



R'^e- 



4r3 



R^ 
2rl 



(93) 



VI. CONCLUSIONS 

In this paper, we have examined the dynamics of ex- 
ponential flights and their relation with the linear Boltz- 
mann equation, a subject that arises in many areas of 
Physics or Biology. In particular, we have focused on 
i) the propagator ^(r|n), which describes the ensemble 
evolution of the transported particles as a function of the 
number of collisions, and ii) the collision density ^(r), 
which is related to the particle equilibrium distribution. 



Moreover, the connection between the number of colli- 
sions and time has been examined. We have provided 
the framework for a generic d-dimensional setup, which 
allows emphasizing the key role of d in determining the 
properties of 5'(r|ri) and ^'(r). 

In this context, we have shown that knowledge of ^'(r) 
formally allows deriving the moments of the residence 
time (or equivalently of the collision number) within a 
given volume, which is key in assessing many physical 
properties of the system at study. The role of bound- 
ary conditions has been explored by considering leakages 
from a given domain, via the method of images. The 
behavior of the spatial moments of the particle ensemble 
has been examined, as well. 

We have then provided specific results for one-speed 
isotropic transport in infinite as well as bounded do- 
mains, and for absorbing or purely scattering media. The 
case d = 1 has been considered as a prototype model of 
exponential flights along a straight line, where only two 
directions (forward or backward) are possible. Due to 
this simplification, most quantities can be explicitly de- 
rived. The case d — 2 has been analyzed in detail: despite 
the calculations being non trivial, in some cases closed- 
form results can be obtained. In particular, we have 
provided an expression for the collision density, which, 
coupled with the method of images, might be useful for 
a realistic description of migration on bounded surfaces. 
The case d = 3 is key in most real-world applications, 
such as the propagation of neutrons or photons in matter, 
or electrons in bulk semiconductors. Unfortunately, this 
case turns out to be hardly amenable to closed-form ana- 
lytical formulas, and most results concern the asymptotic 
behavior of the particles, either close to or far from the 
source. Finally, the case d — 4 has been considered for 
the sake of completeness. Moreover, Monte Carlo simula- 
tions have been performed so as to validate the proposed 
results and support the analysis of the asymptotic be- 
havior. A good agreement is found between theoretical 
predictions and numerical simulations. 

By virtue of the increasing power of Monte Carlo meth- 
ods in solving realistic three-dimensional transport prob- 
lems, one might argue that finding closed- form results for 
simple systems has a limited interest. However, we are 
persuaded that analytical and asymptotic formulas may 
turn out to be useful, in that they can help in improv- 
ing Monte Carlo algorithms by a clever use of the rela- 
tion existing between different variables. Moreover, as 
exponential flights are a transversal field, the cross-over 
between distinct areas of science might hopefully shed 
some light at achieving long-standing goals in transport 
theory, such as full analytical solutions for 3d systems. 
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Appendix A: Spatial moments 



Then, by remarking that for even n 



We begin by computing the m-th coefficient of the 
Taylor expansion of z^^'^^^Jd/2-ii'''z) with respect to z, 
which reads 



1 9™ 
ml dz™ 



(i) 



l+m-l 



z=0 



r(i + f)r(^) 



(Al) 

for even m, and zero otherwise. Apply now the m-th 
derivative to a function f{z) such that f{r) has a spher- 
ical symmetry. Recalling then the definition of the mo- 
ment (r™) from Eq. (|35|) . we have 



(A2) 



1 a™ _ i™2i-^-™(27r)3 



r(i + ^)r(^^) fid 



Rearranging the coefficients, we can finally express the 
spatial moments of f{r) in terms of the m-th z-derivative 
of f{z), namely, 



r(f)r(i±^) d{iz) 



(A3) 



Appendix B: Collision density for d = 2 

In order to find the collision density associated to the 
free propagator, we begin by decomposing the sum over 
the collisions n into even and odd index, i.e., 



*(r) = lim 



\I'(r|n) + ^ *(r|n) 



3dd 



(Bl) 



n even 

we get 



l-{l + z' 



(B2) 



for large N. We have neglected a constant term of the 
kind [log(2) — 7/2]/(27r), which is small compared to 
log(iV). 

For odd n, we can use the series representation 

T / I N ■s-^ 2-^-''^(2-^■)^^-l+^■) ,. 



koven ^r(n-|)r(2-H|-f )r(f) 
2-^-^r(2-t) , 



(B4) 



Now, carrying out the double sum over odd n and over 
k, from Eq. (jB4| we get 



, , log(V7V) + ^+Ez(-r)-log(r) 
*od<iW ~ , (B5) 

for large N. Again, to obtain this result we have ne- 
glected a constant term of the kind [log(2) — 7/2] /(27r), 
which is small compared to log(A^). 

Hence, by summing up we finally obtain 

^^r) ^ logW + ^+E^(-^)-21ogW ^ (gg) 
27r 
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